#####################################################
#####################################################
## Seminář o práci v R Studiu - Část třetí
## 
## Nástroje
#####################################################
#####################################################

# Program:
# 1. Konstrukce intervalů spolehlivosti
# 2. Korelace
# 3. T-test
# 4. Kontingenční tabulky




###################################################
## 1. Konstrukce intervalů spolehlivosti
###################################################

##########################
## 1.1 Analyticky
##########################

# V souboru "weight.Rdata" je možné najít hodnoty váhy u 50 lidí.
# Vypočítejte 95% interval spolehlivosti průměru.

# Nahrajme data:
library(foreign)
load(file= "data/weight.Rdata")

# Podívejme se na data:
weight

summary(weight)

hist(weight)

length(weight)

# Standardní chyba průměru vzorku: sd/sqrt(n).
weight.se <- sd(weight)/sqrt(length(weight))
weight.se

# využití normální distribuce
ci.low <- mean(weight) + qnorm(0.025, mean = 0, sd = 1)*weight.se
ci.up <-  mean(weight) + qnorm(0.975, mean = 0, sd = 1)*weight.se
ci.an <- c(ci.low, ci.up)
ci.an # průměr ostatních vzorků bude uvnitř tohoto intervalu v 95 % případů

# využití t-distribuce
ci.t <- mean(weight) + qt(c(0.025,0.975), df = 49)*weight.se
ci.t



###################################################
##  Cvičení I: Jiná spolehlivost.
#   a) Analyticky odhalte s pomocí t-distribuce intervaly spolehlivosti 
#      na úrovni 90 % a 99 %.






mean(weight) + qt(c(0.05,0.95), df = 49)*weight.se
mean(weight) + qt(c(0.005,0.995), df = 49)*weight.se



##########################
## 1.2 Simulace
##########################

# Nastavíme počet simulací:
nsim <- 1000

# Provedeme 1000 výběrů z normálního rozložení
# s průměrem a standardní chybou jako v datech. Proč?
simweight <- rnorm(n = nsim, mean = mean(weight), sd = weight.se)
simweight

# Zakresleme histogram:
hist(simweight, col="grey80", breaks=20, freq=FALSE,
     main="Výběrová distribuce průměru vzorku (simulační)", xlab = "Vaha")

# Do grafu nyní přidáme:
# linii pro hustotu rozložení
lines(density(simweight)) 

# výpočet empirických kvantilů z distribuce
ci.sim <- quantile(simweight,c(0.025,0.975))

# zakreslení kvantilů
segments(ci.sim[1], 0, ci.sim[1], 0.2, col="red",lwd=2, lty = 2)
segments(ci.sim[2], 0, ci.sim[2], 0.2, col="red",lwd=2, lty = 2)

# zakreslení intervalů spolehlivosti z předchozí analytické kalkulace
segments(ci.an[1], 0, ci.an[1], 0.2, col="green",lwd=2, lty = 2)
segments(ci.an[2], 0, ci.an[2], 0.2, col="green",lwd=2, lty = 2)

ci.sim - ci.an # rozdíl by měl být blízko nuly!



###################################################
## 2. Korelace
###################################################
rm(list=ls())

# Importujme data o výsledcích prezidentských voleb v roce 2018 v obcích.
library(readxl)
prezident <- read_excel("data/prezident2018_kor.xlsx")
# zdroj: úprava na základě dat https://volby.cz/opendata/prez2018/prez2018_opendata.htm 

# Jaký byl vztah mezi ziskem Jiřího Drahoše v 1. a 2. kole?
cor(prezident$prez18_1k_drahos_proc, prezident$prez18_2k_drahos_proc)


###################################################
##  Cvičení II: Vztah zisků Miloše Zemana a Pavla Fischera.
#   a) Vypočítejte korelační koeficient vztahu zisku Pavla Fischera v 1. kole 
#      a Miloše Zemana ve 2. kole.
#   b) Interpretujte toto zjištění.





cor(prezident$prez18_1k_fischer_proc, prezident$prez18_2k_zeman_proc)


# Zobrazit můžeme hned několik korelačních vztahů najednou.
# Nejprve proto vytvoříme korelační matici.
help(cor)
correl <- cor(prezident[,2:12])
install.packages("ggcorrplot")
library(ggcorrplot)
p.mat <- cor_pmat(prezident[,2:12]) # spočítáme také p-hodnoty jednotlivých vztahů
# Chceme totiž odhalit vztahy, které nejsou statisticky signifikantní.

# Nyní můžeme zakreslit všechny výsledky.
ggcorrplot(round(correl,1), hc.order = F, type = "lower", p.mat = p.mat, 
           insig = "blank", lab = F, colors = c("#DC1819", "white" ,"#000A50"), 
           method = "circle", show.diag = F, outline.color = "grey", 
           title = "Prelivy hlasu v prezidentskych volbach")

# Interpretujte zjištěné výsledky



###################################################
## 3. T-test
###################################################

rm(list=ls())

# Soubor fotp.xlsx obsahuje informace o skóre svobody tisku organizace Freedom House.
library(readxl)
data <- read_excel("data/fotp.xlsx", sheet = 2, na = "-")
# zdroj: upraveno na základě https://freedomhouse.org/report/table-country-scores-fotp-2017

# Z datasetu vyřadíme země, kde neproběhlo měření ve všech letech.
data <- data[!is.na(rowSums(data[,2:10])),]

# Proveďte test, zda je rozdíl ve skóre mezi lety 2016 a 2015 statisticky signifikantní.
t.test(data$`2016`, data$`2015`, paired = TRUE, alternative = "two.sided")


# Jelikož je p-hodnota nižší než 0.05 (95% úroveň jistoty), je rozdíl signifikantní.
# Průměrné skóre bylo v roce 2015 o téměř půl bodu nižší než v roce 2016.
# V roce 2016 tak došlo ke zhoršení svobody tisku oproti předchozímu roku.


# Zkuste si změnit parametry funkce t.test - změňte úroveň jistoty a uvidíte, že p-hodnota zůstává stejná,
# ale samozřejmě se mění intervaly spolehlivosti. Sledujte, kdy obsahují nulu a kdy nikoliv 
# (to je totiž přímo spojeno právě s p-hodnotou).



###################################################
##  Cvičení III: Srovnání delší časové řady.
#   a) Máte hypotézu, že v roce 2016 se oproti roku 1996 zhoršila úroveň svobody tisku.
#   b) Jaký je rozdíl mezi oběma roky?
#   c) Jaká je pravděpodobnost platnosti vaší hypotézy?





t.test(data$`2016`, data$`1996`, paired = TRUE, alternative = "greater")
# V roce 2016 bylo průměrné skóre vyšší o 3,53 bodu oproti roku 1996.
# Alternativní hypotéza je platná na 99,97 %.



###################################################
## 4. Kontingenční tabulky
###################################################

rm(list = ls())

# Existuje vztah mezi úrovní vzdělání a důvěrou v prezidenta republiky?
data <- read.csv("data/V1809_F1.csv") # data CVVM ze září roku 2018
# zdroj: Český sociálněvědní datový archiv http://nesstar.soc.cas.cz/webview/ 

# Vytvoření nových předmětných proměnných pro jednoduchost:
data$duvera_prez <- data$PI.1A 
data$vzdelani <- data$t_VZD

# Zmenšení datového setu:
data <- data[,c("duvera_prez","vzdelani")]

# Kontrola podoby proměnných:
table(data$vzdelani)
table(data$duvera_prez)
data <- data[data$duvera_prez!=9,] # vynechání osob, které odpověděli, že neví

# Rekódování proměnných:
data$duvera_prez <- as.factor(data$duvera_prez)
data$vzdelani <- as.factor(data$vzdelani)
library(plyr)
data$duvera_prez <- revalue(data$duvera_prez, c("1"="rozhodne duveruje", "2"="spise duveruje",
                                                "3"="spise neduveruje", "4"="rozhodne neduveruje"))
data$vzdelani <- revalue(data$vzdelani, c("1"="zakladni", "2"="vyuceni",
                                          "3"="stredni", "4"="vysokoskolske"))

# Vytvoření kontingenční tabulky:
tabulka <- table(data$duvera_prez,data$vzdelani)

# Chí-kvadrát test:
chisq <- chisq.test(tabulka)
chisq

# Porovnání pozorovaných a očekávaných hodnot:
chisq$observed  
chisq$expected

# Zakreslení reziduí chí-kvadrát testu jednotlivých buněk:
install.packages("corrplot")
library(corrplot)
corrplot(chisq$residuals, is.cor = FALSE) # pozitivní rezidua jsou modrá, zatímco negativní červená

# Výpočet podílu jednotlivých vztahů na celkovém výsledku chí-kvadrát testu:
podil <- 100*chisq$residuals^2/chisq$statistic
round(podil, 3)
corrplot(podil, is.cor = FALSE)  




